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Oh' We extend existing semi-analytic models of galaxy formation to track atomic and 

q molecular gas in disk galaxies. Simple recipes for processes such as cooling, star for- 

■ mation, supernova feedback, and chemical enrichment of the stars and gas are grafted 
' on to dark matter halo merger trees derived from the Millennium Simulation. Each 

d , galactic disk is represented by a series of concentric rings. We assume that surface 

density profile of infalling gas in a dark matter halo is exponential, with scale radius 
rd that is proportional to the virial radius of the halo times its spin parameter A. As 
^ , the dark matter haloes grow through mergers and accretion, disk galaxies assemble 

■ from the inside out. We include two simple prescriptions for molecular gas forma- 
| tion processes in our models: one is based on the analytic calculations by Krumholz, 

(T) • McKee & Tumlinson (2008), and the other is a prescription where the H2 fraction 

fS| ' is determined by the pressure of the interstellar medium (ISM). Motivated by the 

, observational results of Lcroy et al. (2008), we adopt a star formation law in which 

■ Ssfr oc Sh 2 in the regime where the molecular gas dominates the total gas surface 
I density, and Esfr oc £g as where atomic hydrogen dominates. We then fit these mod- 

t— I , els to the radial surface density profiles of stars, HI and H2 drawn from recent high 

kJ ' resolution surveys of stars and gas in nearby galaxies. We explore how the ratios of 

. £^ I atomic gas, molecular gas and stellar mass vary as a function of global galaxy scale 

^ . parameters, including stellar mass, stellar surface density, and gas surface density. We 

$_i ' elucidate how the trends can be understood in terms of three variables that determine 

the partition of baryons in disks: the mass of the dark matter halo, the spin parameter 
of the halo, and the amount of gas recently accreted from the external environment. 
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1 INTRODUCTION 

Before we can reliably compute how galaxies form stars and 
evolve as a function of cosmic time, we must understand the 
physical processes that regulate the balance between neutral 
and molecular gas in their interstellar media. Only if H2 
forms, will gravitationally unstable clouds cool and collapse 
to high enough densities to trigger star formation in the first 
galaxies. It is also generally believed that star formation 



occurs exclusively in molecular clouds in all galaxies at all 
epochs. 

Many galaxy formation models adopt the so-called 
"Kennicutt-Schmidt" law (hereafter K-S law, Schmidt 1959, 
Kennicutt 1998) to prescribe the rate at which a disk galaxy 
of given cold gas mass and scale radius will form its stars. 
This has the form 

Esfr oc E" as (1) 
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where Esfr represents the star formation rate surface den- 
sity, Ega S is the total surface density of the cold gas in the 
disk, and the exponent n = 1.4 is often adopted. Some 



2 Fu et al. 



semi-analytic models also account for a critical density be- 
low which disks become gravitationally stable and star for- 
mation no longer occurs (e.g Kauffmann 1996, De Lucia & 
Blaizot 2007). In this case, 

E S FR oc [E 

gas S cr it 

] (2) 

where the critical density E cr it is evaluated using the disk 
stability criterion given in Toomre (1964). In both Equa- 
tions HI and <j2j , the star formation rate surface density 
is proportional to the total surface density of cold gas (i.e. 
both HI and H2 components) in the galaxy. This prescrip- 
tion was motivated by the analysis of 97 nearby galaxies 
by Kennicutt (1998), which showed that star formation is 
more tightly correlated with E gas than with Eh 2 ■ There have 
been studies in apparent disagreement with these conclu- 
sions; for example, Wong & Blitz (2002) found that the re- 
lation between Espr and Eh 2 is stronger than that between 
Esfr and E gas in galaxies with high molecular gas frac- 
tions. In recent years, high quality, spatially-resolved maps 
of the cold gas have become available for samples of a few 
dozen nearby galaxies. Examples of such data include HI 
maps from The HI Nearby Galaxy Survey (THINGS) and 
CO maps from the Berkeley-Illinois-Maryland Association 
Survey of Nearby Galaxies (BIMA SONG) and HERA CO- 
Line Extragalactic Survey (HERACLES). Measurements of 
the rate at which stars are forming at different radii in the 
galaxy are provided by Spitzer and GALEX observations. 
The combination of these different data sets has led to im- 
portant new constraints on the relationship between star 
formation and gas in galactic disks. Bigiel et al. (2008) stud- 
ied 18 disk galaxies and showed that H2 forms stars at a 
roughly constant efficiency in spirals at radii where it can 
be detected. Their results suggest a star formation law of 
the form 

Esfr oc E^* 0,2 (3) 

Motivated by these findings, galaxy formation model- 
ers are now progressing beyond a simple single-component 
view of the cold phase of the interstellar medium, and are 
attempting to model the formation of molecular hydro- 
gen in galaxies. Gnedin, Tassis & Kravtsov (2009) included 
a phenomenological model for H2 formation in hydrody- 
namic simulations of disk galaxy formation. Their model 
includes nonequilibrium formation of H2 on dust and ap- 
proximate treatment of both its self-shielding, and shielding 
by dust from the dissociating UV radiation field. Dutton 
(2009) and Dutton & van den Bosch (2009) utilized the 
empirically-motivated hypothesis of Blitz and Rosolowsky 
(2004,2006) that hydrostatic pressure alone determines the 
ratio of atomic to molecular gas averaged over a particular 
radius in the disk in their analytic models of disk forma- 
tion in a ACDM cosmology. They analyzed the radial dis- 
tribution of stars and star formation in their disks, but did 
not focus very much on gas properties in their model. There 
have also been some attempts to predict the balance between 
atomic and molecular gas in galaxies at different redshifts 
by post-processing the publically available outputs of semi- 
analytic galaxy formation models (e.g. Obreschkow et al. 
2009). This work also used the same Blitz and Rosolowsky 
(2004, 2006) prescription to predict the fraction of molecu- 
lar gas in disks. However, the Obreschkow et al. approach 
is not self-consistent, because the simulations have been run 



assuming a "standard" Kennicutt-Schmidt law for star for- 
mation and the presence or absence of molecular gas has 
no influence on the actual evolution of the galaxies in the 
model. 

In this paper, we develop new semi-analytic models 
that follow gas cooling, supernova feedback, the assembly of 
galactic disks, the conversion of atomic gas into molecular 
gas as a function of radius within the disk, and the conver- 
sion of the gas into stars. In the 1990's, semi-analytic mod- 
els of galaxy formation were developed into a useful tech- 
nique for interpreting observational data on galaxy popula- 
tions (e.g. Kauffmann, White & Guiderdoni 1993; Cole et al. 
1994; Somerville & Primack 1999). In the first decade of the 
new Millennium, considerable effort went into grafting these 
models on to large N-body simulations of the dark matter 
component of the Universe. These efforts began with rel- 
atively low resolution simulations (Kauffmann et al. 1999), 
but have rapidly progressed to simulations with high enough 
resolution to follow the detailed assembly histories of mil- 
lions of galaxies with luminosities well below L, (Croton et 
al. 2006; Bower et al. 2006; De Lucia & Bliazot 2007; Guo 
et al. 2010). 

Our new models are an extension of the techniques de- 
scribed in Croton et al. (2006) and De Lucia & Blaizot 
(2007) and are implemented using the merger trees from 
the Millennium Simulation (Springel et al. 2005). We ex- 
plore two different "recipes" for partitioning the cold gas 
into atomic and molecular form: a) a prescription based on 
the analytic models of H2 formation, dissociation and shield- 
ing developed by Krumholz, McKee & Tumlinson (2009), in 
which the molecular fraction is a local function of the sur- 
face density and the metallicity of the cold gas, b) the same 
pressure-based formulation explored by Obreschkow et al. 
(2009). 

We first use our models to calculate the HI, H2, stellar 
mass and SFR surface density profiles of disk galaxies that 
form in dark matter haloes with circular velocities ~ 
200 km/s (i.e. galaxies comparable to the Milky Way) and 
we compare our results to the THINGS/HERACLES obser- 
vations presented in Bigiel et al. (2008). 

We then turn to the issue of the predicted scaling rela- 
tions between atomic gas, molecular gas and stars for an 
ensemble of disk galaxies forming in dark matter haloes 
spanning a range of different circular velocities. We cur- 
rently enjoy a rich and diverse array of scaling laws that 
describe the stellar components of galaxies. For example, 
the Tully-Fisher relation and the size-mass relation for local 
spiral galaxies play a crucial role in constraining current the- 
ories of disk galaxy formation. Likewise, the scaling laws of 
bulge-dominated galaxies (the Fundamental Plane) provide 
important constraints on how these systems may have as- 
sembled through merging. In contrast, few well-established 
scaling laws exist describing how the cold gas is correlated 
with other global physical properties of galaxies. Surveys 
of atomic and molecular gas in well-defined samples of a 
few hundred to a thousand galaxies are currently underway, 
and this paper will explore what can be learned about disk 
galaxy formation from the results. 

Our paper is organized as follows. In section 2, we briefly 
describe the simulation used in our study as well as the 
semi-analytic model used to track the formation of galaxies 
in the simulation. In section 3, we describe the new aspects 
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of the models presented in this paper, including our spa- 
tially resolved treatment of disk formation in radial bins, 
the recipes that prescribe how atomic gas is converted into 
molecular gas, and our new prescriptions for star formation 
and feedback. In section 4, we compare the radial profiles in 
our models to observations from the THINGS/HERACLES 
surveys, and present the global gas properties of the galaxies 
in our model, such as atomic and molecular gas mass func- 
tions. In section 5, we introduce a set of scaling relations for 
the atomic and molecular gas fractions of galaxies and we 
clarify which aspects of the input physics are responsible for 
setting the slope and the scatter of these relations. Finally, 
in section 6 we summarize our work and discuss our findings. 



2 THE SIMULATION AND SEMI-ANALYTIC 
MODEL 

In this section, we give a brief description of the Millen- 
nium Simulation and the physical processes treated in the 
semi-analytic galaxy formation code L-Galaxies. In the next 
section, we describe our own changes to the code, which in- 
clude a resolved model for disk assembly and new recipes to 
treat molecular gas formation, star formation and supernova 
feedback. 



2.1 Mass skeleton: the Millennium simulation 

The Millennium simulation (Springel et al. 2005) is a very 
large N-body cosmological simulation with N=2160 3 « 10 10 
collisionless particles in a comoving box of 500/i _1 Mpc on 
a side. The mass of each particle is 8.6 x 10 s MqH' 1 . The 
cosmogony is ACDM with parameters Qa = 0.75, Q m = 
0.25, Qbaryon = 0.045, cr 8 = 0.9 and h = 0.73. The outputs 
of the Millennium Simulation are stored in a series of 64 
snapshots (0-63); the redshifts of snapshots 4 to 63 are given 
by the expression 

_ 1n (n-63)(n-28)/4200 
Z n — 1U 



1 



4,5-63 (4) 



The redshifts for snapshots to 3 are zq — 127, z\ = 80, 
Z2 — 50, Z3 = 30 respectively. Snapshot 63 corresponds to 
redshift and the time interval between two snapshots is 
approximately 200 Myr. 

In addition to the main Millennium simulation, there 
is a smaller version with the same cosmological parameters 
and mass resolution, the so-called milli-Millennium, which 
includes N=270' ! particles in a comoving box of 62.5/i _1 Mpc 
on a side. The milli-Millennium simulation has the same res- 
olution to the full Millennium simulation and it is very useful 
for fast exploration of parameter space. The full simulation 
is needed to build up sufficient statistics to accurately char- 
acterize the distributions of galaxy properties in a multi- 
dimensional parameter space. 

2.2 The semi-analytic model: L-Galaxies 

L-Galaxies is the semi-analytic code written by Volker 
Springel, described in detail in Croton et al. (2006) and up- 
dated in De Lucia & Blaizot (2007) and subsequent papers 
by the Munich group. It operates on the dark matter halo 
and subhalo merger trees constructed from the 64 snapshots 
of the Millennium Simulation and specifies the treatment 



of the following physical processes: reionization, gas infall 
and cooling, star formation and metal production, supernova 
feedback, galaxy mergers and star bursts, black hole growth 
and AGN feedback. Here we will briefly review those aspects 
of the models most relevant to the analysis presented in this 
paper. 

As described above, the time interval between two suc- 
cessive snapshots from the Millennium Run is about 200 
Myr. In order to model the cooling and star formation ac- 
curately, the time interval between snapshots is divided into 
20 time steps of around 10 Myr, a time interval that is well 
matched to the evolutionary time-scale of massive stars. All 
the physical processes treated by L-Galaxies are computed 
at each time step, but the properties of dark matter haloes 
are only updated at the beginning of each snapshot. 

The gas cooling processes in L-Galaxies follow the treat- 
ment first outlined in White & Frenk (1991). In each dark 
matter halo, the hot gas is assume to be distributed isother- 
mally with a density profile 



TTthot 

4irR vil r 



(5) 



where i? v ir is the virial radius of a halo, defined as the ra- 
dius within which the dark matter density is 200 times the 
critical density, and rrihot is the mass of hot gas within R v i r . 
The local cooling time of hot gas is the ratio of the specific 
thermal energy to the cooling rate per unit volume 



tc 



(r) = 



3[j,m p kBT 
2 Ps (r)A(T, Z) 



(6) 



where fim p is the mean particle density, fee is the Boltz- 
mann constant, and A(T,Z) is the cooling rate. A(T,Z) 
is dependent on the metallicity and the virial temperature 
T = 35.9 (Kir /km s" 1 ) 2 K of the hot halo gas. We adopt 
the cooling function computed by Sutherland & Dopita 
(1993) and the cooling radius r coo i is defined as the ra- 
dius where t coo \ is equal to the dynamical time of the halo 
^dyn — Rvir/vvir = 0. 1H (z) ~ 1 . The cooling rate can be writ- 
ten as 



dm c 



dt 



■ 47rp g (r coo i)r cool 



dr cool 
dt 



As discussed in White & Frenk (1991), when r coo i < 
Rvir, the gas is expected to cool quasistatically from the hot 
gas halo. The mass that cools out at each time step may be 
written as 



Am c 



0.5m hot r cool v vil R vil At, 



(8) 



where At is the time step. If r coo i > R v i T , the halo is in the 
rapid cooling regime, and all the halo gas that has not al- 
ready condensed on to galaxies will accrete on to the central 
object in a so-called "cold flow" . 

Following Mo et al. (1998), the scale length of the disk 
rd that forms through cooling at the centre of a dark matter 
halo can be written 

r d = 7^ r vir (9) 

in which A is the spin parameter of the dark matter halo, 
defined as 



(10) 



where J and E are the angular momentum and energy of 
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the dark matter component of the halo. The disk dynamical 
time is given by 

fdyn = 3r d /Uvir- (11) 

In galaxy disks, cold gas is converted into stars using 
the formula 

JO (?n gas < merit) . . 

m* = < (12) 
la (m gas — m cr it) /tdyn (m gas ^ merit) 

in which the constant a is the star formation efficiency, and 
m C rit is the critical mass of gas in the disk. Following Kenni- 
cutt (1989) and Kauffmann (1996), the critical density can 
be approximately expressed as 

E-» = 0J»^Pc-"(^b r )(^)" 1 (13) 

and the critical mass m cr it (assuming an outer disk radius 
of 3rd and a fiat gas density profile) is therefore 

merit = 5.7 x 1O 6 M (^) (g-) . (14) 

We note that when m gaa < m cr it, star formation stops over 
the entire disk. The reason why the star formation rate is 
taken to be proportional to m gas — m cr it when m gas > m cr it, 
is to ensure that the star formation rate is a continuous 
function of the gas mass. 

Energy from supernova explosions can reheat the cold 
gas in disks to hot gas, and it can also eject some of the hot 
gas from the halo. The amount of reheated cold gas is 

Am re h ea t = EdiskAm,, (15) 

where Am, is the mass of stars formed in a given time step, 
and edisk is the supernova reheating efficiency. The total en- 
ergy released into feedback processes is approximately 

ASsn = 0.5e ha ioAm,Dg N , (16) 

where 0.5vg N is the energy in supernova ejecta per unit 
mass of newly formed stars, and eh a io is the efficiency with 
which the supernovae convert disk gas to halo hot gas. 
In Equation (|16[) . «sn = 630km s -1 is adopted, based on 
standard supernova theory. If the excess energy in hot gas 
ASexcess = AiJsN — Ai5rche a t > 0, part of the hot gas will 
be ejected from the halo. The amount of ejected hot gas is 

2A_E eX cess A ^SN a f-tT\ 

Am e ject = 2 = £ haioAm.- Am rc hc a t (17) 

^'vir ^ 7 vir 

The ejected gas will not be lost permanently. As the 
dark matter halo grows through accretion and merging, the 
ejected gas will be reincorporated into the central halo and 
become part of the hot phase once again. At each time step 
(At), the amount of reincorporated gas is 

Am r eincorpor a te — 7ej (^ejeotedAdyn ) At (18) 

where m e j C cted is the mass of gas currently in the ejecta 
reservoir and 7 e j is the fraction of the ejected gas that is 
reincorporated per halo dynamical time. The default value 
of 7 c j is 0.5, meaning that half the ejected gas will be rein- 
corporated into the central halo after one dynamical time 
(where = Rvir/v vi r) 

The Instantaneous Recycling Approximation (IRA) is 
adopted to model the chemical enrichment in the galaxies. 



Assuming the initial mass function of Chabrier (2003), the 
yield is 0.03 and the fraction of the stellar mass returned 
to the cold interstellar medium is 43%. The flow of metals 
through the cold, hot and ejected components is tracked in 
exactly the same way as the flow of gas. 

In the L-Galaxies code, galaxy mergers are divided into 
two classes : minor mergers and major mergers. The bary- 
onic mass (gas+star) ratio of the satellite galaxy to the 
central galaxy m sa t/m C en determines the merger type. In 
a minor merger (m sa t/m C en < 0.3), the stars of the satellite 
galaxy are added to the bulge of the central galaxy, while 
the cold gas of the satellite galaxy is added to the disk of 
central galaxy. In a major merger (m sa t/m C en > 0.3), the 
stellar disks of both galaxies are destroyed and a bulge is 
formed. In both minor and major mergers, the mass of stars 
formed in the associated starburst is (following Somerville 
et al. 2001) 

Am, = /3 b urst (m sat /m ccn ) Qburst m gas . (19) 

In minor mergers, the stars formed in the burst are added to 
the disk, but in major mergers all stars formed in the burst 
are added to the bulge. 

Detailed descriptions of other processes such as reion- 
ization, black hole growth and AGN feedback can be found 
in Croton et al. (2006), and we inherit these recipes without 
modification. 



3 MODIFICATIONS TO THE MODELS 

If we wish to study the relations between neutral gas, molec- 
ular gas and star formation in galaxies, a single zone treat- 
ment of the gas and stars is not sufficient. In this section, 
we will describe how we have changed the code so that it is 
able to track the radial distribution of gas and stars in galac- 
tic disks. Once we are able to specify the surface density of 
stars and gas as a function of radius in the disk, it is then 
a simple matter to incorporate the prescriptions of Blitz & 
Rosolowsky (2006) or Krumholz et al. (2009) to predict the 
fraction of the gas that is in molecular form. 

3.1 Modelling the radial profiles of the disks 

In the standard version of the L-Galaxies code, we track the 
total stellar mass and gas mass as a function of cosmic time 
for each galaxy in the simulation. In our new implementa- 
tion, we adopt a fixed set of 30 radial "rings" for each galaxy, 
and follow the build-up of stars and gas within each ring. 
The radius of each ring is given by the geometric series 

n = 0.5 x 1.2* [/i _1 kpc] (i = 1, 2.. .30) (20) 

where innermost ring has a radius of 0.8 kpc and the outer- 
most ring has radius ~ 160 kpc. We have experimented with 
other ways of setting up the rings, for example, arithmetic 
radial division, and with changing the number of rings. The 
radial profiles are insensitive to the precise scheme, so long 
as the adopted number of rings is sufficiently large. 

Following Mo et al. (1998), we assume the gas that cools 
at a given time step is distributed exponentially with surface 
density profile 

S gas (r) = E gas exp (-r/nnfali) (21) 
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The radial distribution of The cumulative radial The radial distribution of the gas, including 

the gas accreted at each output time distribution of the accreted gas star formation and supernova feedback 




r [kpc] r [kpc] r [kpc] 

Figure 1. Illustration of how our superposition method works. Results are shown for a "Milky Way" type galaxy in a halo with present- 
day circular velocity v viI ~ 200 km/s . In the left panel, we show the radial profiles of the infalling gas for the "main progenitor" of the 
galaxy at each output time in the simulation. Each profile is colour-coded according to the output time, as indicated in the colour bar 
on the right. In the middle panel, we plot the cumulative profile of all the gas that has been accreted up to each output time. In the 
right panel, we show the actual gas profile of the galaxy at each output time, after star formation and supernova feedback processes have 
been taken into account. 



If we assume that the angular momentum of the infalling gas 
is conserved, the scale length of the infalling gas is ri n f a n = 
r c i = (A/V2) r v i r , where r v j r is the virial radius of the halo 
and r<j and A have the same meaning as in Equations ([9| & 
(|10p . In Equation (|21[1 . Eg as is the central surface density of 
the infalling gas and is given by 

^ gas _ (22 > 

where m CO oi is the mass of gas that cools in a given time 
step. Note that we assume that metals that leave the galaxy 
disk are uniformly mixed with the halo gas, which results in 
a uniform metallicity for the infalling gas. 

One important issue that must be addressed is how the 
gas that accretes at some time t% is connected to the gas that 
accretes at a later time ti+i. We adopt the simplest possible 
superposition scheme, whereby the pre-existing radial profile 
of the gas remains unchanged and the new infalling gas is 
superposed directly on to it. 

An illustration of this scheme is presented in Fig. [T] The 
figure shows the gas accretion history of a "Milky Way" type 
disk galaxy with stellar mass of M* = 1O 1O ' 7 M0 at z=0, re- 
siding at the centre of a dark matter halo with virial velocity 
tVir = 207km/s. The left panel shows the radial profiles of 
gas that accretes on to the main progenitor of the galaxy 
at each snapshot. The middle panel shows the cumulative 
profile of all the gas that has been accreted up to that snap- 
shot. The right panel shows the actual gas profiles at each 
snapshot, after star formation and supernova feedback are 
included. Each curve is colour-coded according to the time 
elapsed since the Big Bang: red represents early times and 
blue represents late times. 

At early times, the disk is small and compact. As the 
Universe evolves, the disks grow in size. The accretion rate of 
gas on to the disk reaches a maximum around 4-8 Gyrs after 
the Big Bang and then declines until the present day. Q In 

1 This is similar to the evolution of the global star formation 



the right panel, we can see how the gas in the central region 
of the disk is depleted at late times as the gas is transformed 
into stars. We note that similar "inside-out" disk formation 
models have been explored in the past (e.g. Kauffmann 1996; 
Dalcanton, Spergel & Summers 1997; Van den Bosch 1998; 
Avila-Reese, Firmani & Herandez 1998; Boissier & Prantzos 
1999; Fu et al. 2009). 

When two galaxies merge, we simply add the stellar 
and gas radial profiles of the satellite galaxy to those of the 
central galaxy and then allow for a star burst in each radial 
bin in the same way as in the original L-Galaxies code. This 
is probably not an accurate description of what happens 
in reality, but our intention in this paper is to concentrate 
on the predicted gas profiles and scaling relations of disk- 
dominated galaxies, which do not experience many mergers. 

In this paper, we will not consider the radial distribu- 
tion of the dark matter. We adopt an approximation of a 
constant circular velocity for the whole disk « c ir = v V ir when 
comparing the model results to the observations. 



3.2 The conversion of atomic to molecular gas 

Having specified how the infalling gas is distributed across 
the disk, we are now ready to consider molecular hydrogen 
formation processes. 

We begin by noting that the total cold gas content of the 
galaxy predicted by our models includes both hydrogen and 
helium. The correction factor between cold hydrogen mass 
and total cold gas mass is a factor of 1.4 (Arnett 1996). We 
note that a correction for helium is applied in some obser- 
vational analyses, but not in others, so we must be careful 
to make sure we always compare the models with the data 



rate density of the Universe (Madau, Pozzetti & Dickinson 1997), 
which reaches a maximum around redshift 0.5-1.5, and then de- 
creases. 
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in a self-consistent way. We also note that recent observa- 
tions (e.g Reynolds 2004) of the Local Interstellar Clouds 
in the Milk Way disk show that about 1/3 of all the gas in 
the Milky Way disk is in a warm, ionized phase. Following 
Obreschkow et al. (2009), we include a warm phase correc- 
tion into our models. The warm phase coefficient is £ = 1.3. 
If there is no note to the contrary, the gas surface densities 
are always divided by £ in this paper. 

We have implemented two different molecular gas for- 
mation prescriptions in our models. The first prescription is 
from Krumholz et al. (2009) (hereafter K09), in which the 
H2 fraction is a function of local cold gas surface density and 
metallicity (hereafter H2 prescription 1). The other prescrip- 
tion has its origin in papers by Elmegreen (1989 & 1993), 
Blitz & Rosolowsky (2006), and Obreschkow et al. (2009), 
which have all suggested that the H2 fraction is a function 
of the pressure of the ISM (hereafter H2 prescription 2). 



3.2.1 The H2 prescription of Krumholz et al. 

Krumholz et al. (2008) considered the location of the atomic- 
to-molecular transition in a uniform spherical gas cloud 
bathed in a uniform, isotropic, dissociating radiation field. 
The main result of their calculation is that the amount of 
atomic material required to shield a molecular cloud against 
dissociation by the interstellar radiation field is character- 
ized by two parameters: tr, a measure of the dust opti- 
cal depth of the cloud, and Xi a normalized radiation field 
strength, which is defined as the ratio of the rate at which 
Lyman- Werner photons are absorbed by dust grains to the 
rate at which they are absorbed by hydrogen molecules in a 
parcel of predominantly atomic gas in dissociation equilib- 
rium in free space. 

In Krumholz et al. (2009) and McKee & Krumholz 
(2010), this idealized model was applied to atomic-molecular 
complexes in galaxies in an attempt to elucidate the physi- 
cal processes and parameters that determine the molecular 
content in these systems. The paper demonstrates that be- 
cause of the way the density in the cold phase of the atomic 
ISM varies with the interstellar radiation field, \ ~ 1 in 
all galaxies with only a weak dependence on metallicity. 
The existence of a "characteristic" normalized radiation field 
strength, leads to a simple analytic approximation for the 
fraction of mass in an atomic-molecular complex that is in 
the molecular phase solely in terms of the column density 
of the complex and the metallicity of the gas. Krumholz et 
al. (2009) show that the predictions of their models agree 
well with observations of the atomic and molecular content 
of the cold gas in the Milky Way and in nearby galaxies, 
particularly in regions where E gas > lOM0pc~ 2 . 

In our model, we adopt the molecular fraction values 
/H 2 (E gas , [Z/H]g a s) from K09 to calculate the mass of neu- 
tral and molecular gas in each ring, so that 



Eh 2 = /H 2 S gas 
Shi = (1 - /h 2 ) E g 



(23) 



Since the gas surface density defined in K09 is a local 
surface density, while E gas in our models is an azimuthally- 
averaged surface density, we introduce a clumping factor c t 
to account for the fact that the gas in disks is usually orga- 
nized into higher density clumps in spiral structures. In this 



case, /h 2 should be written as 

/h 2 = /H 2 (c f £ gas , [Z/H] gas ). 



(24) 



In our models, cj is left as a free parameter (see section 3.5 
for a description for how the free parameters are set). 



3.2.2 The pressure related H2 prescription 

Blitz & Rosolowsky (2004 & 2006) proposed that pressure 
alone determines the ratio of atomic to molecular gas av- 
eraged over a particular radius in disk galaxies. In these 
papers, the molecular ratio of a galaxy disk was expressed 



7? mol (r) ee M H2 (r) /Mm (r) = [P (r) /P ] ap , 



(25) 



where Po and ap are constants fit from the observations. 

We follow a similar procedure that outlined in 
Obreschkow et al. (2009) to calculate the pressure in our 
model disk galaxies. According to Elmegreen (1989 & 1993), 
the midplane-pressure of the ISM in disk galaxies can be ex- 
pressed as 

P (r) = jGE gas (r) [E gas (r) + /„ (r) E. (r)] , (26) 

where G is the gravitational constant, r is the radius from 
the galaxy centre, and f a (r) is the ratio of the vertical ve- 
locity dispersions of the gas and the stars: 



(r)/<r» (r) 



(27) 



Observations indicate that <r gas is approximately a con- 
stant across the whole disk (e.g. Boulanger & Viallefond 
1992; Leroy et al. 2008). o-*(r) decreases exponentially with 
a scale length twice that of the stellar disk (e.g Bottema 
1993), i.e 



o~*(r) = er* exp (— r/2r«) , 



(28) 



where al is the stellar velocity dispersion at the centre of 
the disk. Substituting the exponential stellar disk profile 
E* (r) = Ejexp(— r/r») and Equation (f28)l into Equation 
(pT)) , we get 



f , \ cr ga s / E2 



fa 



eTM 



(29) 



where /° is the value of / CT (r) at the centre of the disk. The 
mean value of f a (r) for the whole disk is 



f 27rrE„ (r) f a (r) dr _ E2/°(2r«) 



/ 2-n-rE* (r) dr 



E°r| 



4/° (30) 



Adopting the value quoted by Elmegreen (1993), f a ss 0.4, 
we get fa(r) ~ 0.1i/E*/E*(r). Substituting into Equation 
the pressure of ISM may be written 



P(r) = -GE gas (r) 



; (r)+0.iyTU (r) E° 



(31) 



Equation (|31[) can then be substituted into Equation 
(|25[) to obtain the molecular ratio. We adopt Po = 5.93 x 
10" 13 Pa and a P = 0.92 from Blitz & Rosolowsky (2006), 
which were the mean values obtained from fitting Equation 
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(|25[) to a small sample of nearby galaxies. Our final expres- 
sion for the molecular fraction is 



H 2 prescription 1 



H 2 prescription 2 



/h 2 (r) 



1.38 x 10" 



E* M (r) + 0.1E gM (rWE.(r)£° 



(32) 



where E gas (r), E*(r) are the radial gas and stellar surface 
densities, E°, is the central stellar surface density (units 
are Mopc~ ). Thus, in the pressure-based prescription, the 
molecular ratio depends primarily on the gas surface density, 
but the stellar surface density can be important, particularly 
in gas-poor disk galaxies at low redshifts. 

Finally we note that when we apply our pressure re- 
lated prescription to the disks in our simulation, we do not 
apply any clumping factor or warm phase corrections, be- 
cause Equation (|25[1 refers to the total pressure averaged 
over a particular galactocentric radius. 



3.2.3 Comparison between the two H2 fraction 
prescriptions 

Fig. [2] illustrates the main differences between the two H2 
prescriptions described above. On the left, we plot /h 2 vs. 
E gas for the Krumholz prescription; each curve corresponds 
to a different value of the gas-phase metallicity. On the right, 
we plot /h 2 vs. Ega S for the pressure-based prescription; each 
curve corresponds to a different stellar surface density (note 
that we assume a central stellar surface density of E° = 
2OOOAf0pc -2 in Equation (|32[) . which is typical for a spiral 
galaxy like the Milky Way). 

We see that the molecular fraction decreases at low gas 
surface densities in both panels, but the thresholding effect 
is much stronger for the Krumholz prescription than for the 
pressure-based prescription. In addition, the Krumholz pre- 
scription appears to open up the possibility that molecular 
gas formation is strongly suppressed in low mass galaxies, 
which have both low densities and low metallicities. We will 
come back to these points in future work. 



3.3 Star formation 

In a recent paper, Leroy et al. (2008) measured the local 
star formation efficiency (SFE) (i.e the star formation rate 
(SFR) per unit mass in gas) for a small sample of nearby 
galaxies and compared it with expectations from a num- 
ber of proposed star formation laws. Their basic result is 
that in the inner regions of spirals, where the molecular 
gas constitutes a significant fraction of the total cold gas 
content of the galaxy, the SFE of H2 is nearly constant at 
(5.25 ± 2.5) x 10 -10 yr -1 (corresponding to a H2 depletion 
time about 2 x 10 9 yr). Leroy et al. tested whether there 
were additional dependencies on variables such as free fall 
and orbital time-scales, midplane gas pressure, stability of 
the gas disk to collapse, and the ability of perturbations to 
grow despite shear, but they found no effect. 

These results lead us to implement the following star 
formation "law" in our models: 



Esfr = oEh 2 , (33) 
where the molecular star efficiency a is a constant for all 
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Figure 3. The radial surface density profiles of stars, cold gas, 
HI and H2 for galaxies with masses similar to that of the Milky 
Way at z=0. The light blue curves with error bars are taken from 
the data presented in Leroy et al. (2008) and include galaxies 
with circular velocities in the range 200km/s< ti v i r <235 km/s 
and ulge /M* 15%. The red solid and dotted curves are the 
mean and median values from the models, and the black dashed 
curves show the ±1<t deviations about the mean. The panels on 
the left show the model results when H2 prescription 1 is used, 
and the panels on the right are for H2 prescription 2. 



galaxies, and we adopt a — 4.9 x 10 
value in Leroy et al. (2008). 



yr 



similar to the 



3.4 Constraints from the radial profiles 

In Fig. [3] we plot the radial surface density profiles of stars, 
HI, H2 and total cold gas for "Milky Way" sized disk galaxies 
residing in dark matter haloes with circular velocities in the 
range 200-235km/s at redshift 0. We show results for both 
H2 fraction recipes and compare with data from Leroy et al. 
(2008). We find that we are unable to get a good fit to the 
data, even if we adjust all available model parameters. Al- 
though we can get stellar surface density profiles that match 
the observations very well, the HI and total gas surface den- 
sity profiles are too flat. In addition, the H2 surface density 
is more than an order of magnitude too low. The main rea- 
son for this discrepancy is that a star formation law of the 
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H 2 prescription 1 H 2 prescription 2 




Z gas [M Q pc- 2 } £ gas [A/ pc- 2 ] 

Figure 2. Comparison of the results for the two H2 prescriptions. The left panel shows the relation between H2 fractions and local gas 
surface density predicted by the Krumholz et al. models for four different gas mctallicities. The right panel shows the relation between H2 
fraction and local gas surface density predicted by the pressure-based prescription for four different values of the stellar surface density. 
A central stellar surface density = 2000A/qpc — 2 is assumed in the right panel. 



form given in equation ()33[) leads to a gas consumption rate 
that is steeper than the total gas profile, resulting in rapid 
depletion of the gas in the inner part of the disk. 

Recall that the amount of gas in the disk is regulated 
not only by star formation, but also by supernova feedback 
processes, which act to heat the cold gas. The supernova 
reheating rate is proportional to the mass of newly formed 
stars (Equation ll5|) . so Equation (|33|) then implies that the 
net gas consumption rate will proportional to the surface 
density of H2. Because the H2 fraction is a strongly increas- 
ing function of gas surface density, the H2 surface density 
profiles are always steeper than the total gas surface density 
profiles. As a result, the gas consumption rate has a steeper 
dependence on radius than the gas profile itself, and this 
causes the gas profiles to flatten with time. 

One possible solution to this problem is that the gas 
does not remain fixed at a given radius, but flows inwards, 
thus replenishing the gas that is consumed at the centre of 
the galaxy by star formation. We do not consider radial gas 
flows in this paper. Instead, we make two plausible changes 
to the star formation "law" and feedback prescriptions that 
act to bring the gas profiles in better agreement with obser- 
vations. 

First, we are computing /h 2 based on azimuthally- 
averaged column densities, but this approximation becomes 
increasingly poor in the outer parts of galaxies where the 
mean column density and mean molecular fraction are very 
low. In such regions the star formation tends to be domi- 
nated by small isolated regions where the local column den- 
sity is much higher than the azimuthally-averaged value, and 
as a result the molecular fraction is higher. Since we cannot 
easily capture this effect with our azimuthally-averaged pre- 
scription, we instead modify the star formation law based on 
an empirical fit to the behavior of the star formation rate 
surface density in the outer regions of galactic disks (Bigiel 



H 2 prescription 1 



H 2 prescription 2 




~ -1 1 2 -10 12 3 

logic s hi+h 2 [M pc~ 2 ] log 10 S H i+h 2 [M pc" 2 ] 

Figure 4. Black points show the star formation rate surface den- 
sity as a function of the HI+H2 surface density within each radial 
ring for a set of randomly selected galaxies from our simulation 
at z=0. Blue points are observational data taken from Figure 10 
of Bigiel et al. (2008). Results are shown for both H2 fraction 
prescriptions. 



et al 2008). We adopt 



EsFR = 



cv£h 2 (/h 2 ^ 0.5) 
a'E 2 as (/ H2 < 0.5), 



(34) 



where a has the same value as in Equation (|33f) . and a' = 
0.5cv/E gas |/ H =0.5, so that the star formation rate changes 
continuously at the radius where /h 2 = 0.5. In other words, 
our proposed "hybrid" star formation law results in a star 
formation rate that is proportional to the Eh 2 surface den- 
sity in the regions of the disk where H2 is dominant, but 
a Schmidt-Kennicutt type star formation law Esfr oc E™ as 
with n = 2 is adopted in the HI dominant regions. 

In Fig. [4] we plot Esfr versus Ehi+h 2 for a set of ran- 
domly selected galaxies from our simulation at z=0. The 
black points represent gas surface densities and star forma- 
tion rates evaluated in the 30 radial rings in each of the 
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selected galaxies. Results are shown for both H2 fraction pre- 
scriptions. In order to be consistent with the data in Bigiel 
et al. (2008), the helium component is not included in Fig. 
[4] For comparison, we plot the observed relation between 
Espr and Ehi+h 2 from Bigiel et al. (2008). As can be seen, 
our "hybrid" SFR law provides an excellent match to the 
real data over the entire range is total gas surface density 
probed by the observations. A change in slope occurs at 
Shi+h 2 ~ IOMqpc -2 , which marks the division between 
the H2 and HI dominant regions of the galaxy. 

The change in star formation law helps to steepen the 
gas profiles in the outer disk, but does not solve the problem 
of the over-consumption of gas in the central regions of the 
galaxy. In our standard supernova feedback recipe, the mass 
of cold gas reheated by supernovae is proportional to the 
mass of newly formed stars. We might hypothesize that the 
dissipation of energy input by supernovae is more efficient in 
denser regions (Equation [6]), so that less gas is reheated in 
the inner regions of the disk per unit mass of stars that are 
formed. We have chosen to make a correction to Equation 
(|15p so that the reheated mass is inversely proportional to 
the surface density of gas, 

Am re h ea t = Qiisk Am, , (35) 

^gas 

where the coefficient Eo gas is the same for all galaxies. 

In the remainder of this paper, we will always adopt 
Equations (|34[) & (|35[) to describe star formation and super- 
nova reheating, respectively. 

3.5 Model parameters 

In Table [1] we list the values of the new model parameters 
that we introduced: c f , £, P , a P , a, e disk , E 0gas and k A gn- 
The values of the other parameters in the L-Galaxies code 
are unchanged and they can be found in Table 1 in Cro- 
ton et al. (2006). Note that the values of £, Pq and ap are 
adopted directly from observations and are not tunable pa- 
rameters. Cf, a, tdisk, Eogas an d k agn are free parameters in 
our models, which we adjust to fit both the observed stel- 
lar and gas density profiles and the global stellar and gas 
mass functions (see Sec. 4.2). We will now clarify how each 
of these parameters affects our main results: 

cf. As discussed in Sec. I3.2.T1 the clumping factor cor- 
rects for the difference between the local gas surface den- 
sity that is relevant to the Krumholz et al. H2 fraction pre- 
scription, and the azimuthally-averaged gas surface densities 
predicted by the models. Higher values of Cf lead to higher 
predicted H2 fractions. The left panel of Fig. [2] shows that 
a change in a has the greatest effect in the low gas surface 
density regions of the disk. In practice, we tuned Cf so as to 
best reproduce the observed H2 and HI mass functions at 
z = 0. We found that Cf = 1.5 allowed us to fit both the gas 
surface density profiles of Milky Way-type galaxies and the 
gas mass functions at z=0. 

a: Leroy et al. (2008) find an H2 star formation effi- 
ciency a — (5.25±2.5) x 10 -10 yr _1 from their observational 
data. In our models, the value of a controls the amplitude of 
the stellar mass profiles and the stellar mass function. Our 
adopted value of a — 4.9 x 10 -10 yr _1 allows us to fit the 
stellar mass function at z = and is in remarkably good 
agreement with the value found by Leroy et al (2008). 



edisk and Eo g a S : From Equation (|35[) , we see that 
ediskEogas controls the mass of cold gas reheated by super- 
nova explosion. Higher values of edisk Eo g as lead to lower gas 
surface densities and a lower amplitude of the gas mass func- 
tions. We have chosen to fix edisk at a value of 3.5 (the same 
as in Croton et al. 2006). We tune Eo ga s to fit the gas surface 
density profiles and gas mass functions. 

kagn: is the quiescent hot gas black hole accretion rate, 
which controls the fraction of hot gas in massive dark matter 
haloes that condense on to central galaxy that host super 
massive black holes. Higher values of kagn decrease the ac- 
cretion rate of gas from the surrounding halo and suppress 
the number of very high mass galaxies. Because we changed 
the way star formation is treated in the model, we retuned 
kagn to obtain good fit to the stellar mass function at the 
high mass end. 



4 RESULTS 

4.1 Radial surface density profiles at redshift=0 

In Fig. [5] we compare the radial surface density profiles 
of stars, HI, H2 and the total cold gas in disk galaxies 
at redshift z=0 with the observational results of Leroy et 
al. (2008). In the right panels, we show results for galax- 
ies with circular velocities similar to the Milky Way (200 
km/s< v c i r <235 km/s), and in the left panel we show re- 
sults for disk galaxies with somewhat lower circular veloci- 
ties (170 km/s< Heir <200 km/s). As mentioned in Sec. 13.11 
we assume constant disk radial velocity profile u c i r = «vir for 
simplicity. The Leroy et al. (2008) data compilation consists 
of E» profiles derived from data from the Spitzer Infrared 
Nearby Galaxies Survey (SINGS) (Kennicutt et al. 2003), 
Ehi profiles from The HI Nearby Galaxy Survey (THINGS) 
(Walter et al. 2003), and Eh 2 profiles from the HERACLES 
(Leroy et al. 2008) and BIMA SONG (Heifer et al. 2003) 
surveys. Four of the galaxies in the Leroy et al. (2008) sam- 
ple have circular velocities in the range 200 km/s< v clr <235 
km/s: NGC 0628, NGC 3184, NGC 5194, and NGC 3521; an- 
other four have 170km/s< v cil <200km/s: NGC 3351, NGC 
6946, NGC 3627, and NGC 5055. All 8 galaxies are spi- 
ral galaxies with morphological type later than Sb (Hubble 
type index T > 3). Laurikainen et al. (2007) and Weinzirl 
et al. (2009) have quantified the average relation between 
Hubble type and the ratio of the luminosities of bulges to 
galaxy disks. The fraction of the mass of the galaxy in the 
bulge (B/T = Ml Jul£C /M,) for spiral galaxies later than Sb 
is typically lower than 0.15. Thus, we select the model galax- 
ies in dark matter haloes with virial velocities in the range 
200 km/s< w vir <235 km/s and 170 km/s< v vil <200 km/s 
with B/T ^ 15% to compare to the observations. In each 
sub-panel of Fig. [5] the mean profiles for the models that as- 
sume the Krumholz H2 fraction prescription are plotted as 
red solid curves, the median profiles as a red dotted curves, 
and the ±l(j deviations about the mean as black dashed 
curves. The mean profiles for the models that assume the 
pressure-based H2 fraction prescription are plotted as green 
solid curves. 

In Fig. [5] we see that the mean surface density profiles 
of the gas and the molecular gas in the inner regions of the 
disk are reasonably well fit by an "exponential law". This 
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Table 1. The model parameters 



Parameter Description Value 



Cf clumping factor for H2 prescription 1 (Eg, 12 It 1.5 

§ warm phase correction factor (Sec, 13.21 1.3 

Pq the constant of the relation between molecular ratio and ISM pressure (Eg. I25H 5.93 X 10 — 13 Pa 

cep the index of the relation between molecular ratio and ISM pressure (Eq. I25H 0.92 

a star formation efficiency in (Eg. l33l fel34t 4.9 X 10 _10 yr _1 

e disk the supernova reheating rate (Eg, I15l fe l35l) 3.5 

Sogas constant surface density in the supernova reheating recipe (Eg. 1351) 2.1M0pc -2 

kagn Quiescent hot gas black hole accretion rate 2.0 X 10~ 6 MQyr~ 



170km/s < v vir < 200km/s 200km/s < v vir < 235km/s 




Figure 5. The radial surface density profiles of stars, HI, H2, 
and total cold gas for disk galaxies in our models at redshift 0. 
The light blue curves with error bars are from the observational 
data of Leroy et al. (2008). The red curve shows the mean profile 
for H2 prescription 1, while the dashed black curves show the 
±lcr deviations about the mean. The green curves show the mean 
profiles for H2 prescription 2. All the gas surface density profiles 
include the contribution from helium and a correction for the 
warm phase. The right column shows results for disk galaxies 
with rotation velocities similar to the Milky Way (200 km/s< 
■u v ir <235 km/s), and the left column for lower mass disk galaxies 
(170 km/s< j; vir <200 km/s). 



is mostly a consequence of the fact that the infalling gas is 
assumed to have an exponential profile at each time step. 
Fig. [S] shows that the problem with the overconsumption of 
molecular gas in the inner disk discussed in Section \3. 41 has 
been solved by the new implementation of supernova feed- 
back discussed in the previous section. In agreement with 
observations, the molecular gas dominates in the inner disk, 
the neutral gas dominates in the outer disk, and the H2 pro- 
files closely track the stellar profiles. Comparing the red and 
green curves, we see that the two H2 fraction prescriptions 
give very similar results. Note that we have tuned the avail- 
able free parameters to reproduce the observational data as 
well as possible, so this is not entirely surprising. The main 
differences between the two prescriptions occur in the outer 
disks, where Eh 2 falls off less steeply for the pressure-based 
prescription, particularly for galaxies with lower circular ve- 
locities, which have lower metallicities. The results for inner 
disks are almost identical for the two prescriptions. 



4.2 Stellar and gas mass functions at z=0 

In Fig. [S] we show the HI, H2 and the stellar mass functions 
predicted by our models at z = 0. The results are compared 
to the HI mass function derived from the HI Parkes All Sky 
Survey (HIPASS) by Zwaan et al. (2005), the stellar mass 
function derived from the data release 7 of the SDSS by Li 
et al. (2009), and the H2 mass function derived from the 
Five College Radio Astronomy Observatory (FCRAO) Ex- 
tragalactic CO Survey by Keres et al. (2003) assuming a con- 
stant CO-H2 conversion factor (blue), and by Obreschkow 
& Rawlings (2009) assuming a variable CO-H2 conversion 
factor that scales with the B-band magnitude of the galaxy, 
as proposed by Boselli, Lequeux & Gavazzi (2002) (black). 
The red curves in Fig. |S] are our model results with H2 pre- 
scription 1, while the green curves show the results for H2 
prescription 2. 

Note that we only plot the stellar mass function down to 
a limiting mass of M» = 1O 9 ' 5 M0, for which the Millennium 
Simulation has sufficient resolution to accurately track the 
formation history of the host halo. The stellar mass func- 
tions predicted by the models agree very well with the data 
for galaxies with stellar masses greater than 1O 1O M0. The 
models produce too many galaxies with masses less than 
10 10 Mq. Guo et al. (2010) have implemented the galaxy 
formation model of De Lucia & Blaizot (2007) in a higher 
resolution simulation, and show that the mismatch with the 
faint end of the galaxy mass function becomes increasingly 
severe as one goes to lower masses. Since we are not con- 
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Figure 6. The mass functions of HI, H2 and stars for galaxies at redshift zero. The red and green curves are model results for H2 
prescription 1 and 2, respectively. The HI and H2 mass functions include the correction for the warm phase but do not include the helium 
component. The open symbols with error bars indicate the observational data, while the blue and black curves indicate the best-fit 
Schechtcr functions. The stellar mass function is taken from Li et al. (2009) analysis of the final data from the Sloan Digital Sky Survey. 
The HI mass function was derived from HIPASS by Zwaan et al. (2005). The H2 mass function is from the FCRAO Extragalactic CO 
Survey; the blue circles and curves assume a fixed CO-H2 conversion factor (Keres, Yun & Young 2003), while the black curves assume 
a variable CO-H2 conversion factor as given in Obrcschkow & Rawlings (2009). 



cerned with dwarf galaxies in this paper, we will ignore this 
problem for the moment. 

The agreement with the HI mass function of Zwaan et 
al. (2005) is excellent. It is also interesting that the H2 mass 
function predicted by the model is in closer agreement with 
the observational results derived using a constant conversion 
factor. 

Finally, we note that the stellar and cold gas mass func- 
tions predicted by the models at z=0 are insensitive to the 
choice of H2 fraction prescriptions. 



5 GLOBAL SCALING RELATIONS BETWEEN 
THE MOLECULAR GAS, NEUTRAL GAS 
AND STELLAR MASSES OF GALAXIES 

In this section, we investigate the correlations between the 
atomic gas fraction Mhi/M*, the molecular gas fraction 
Mh 2 /M, and the molecular-to-atomic gas ratio Mh 2 /Mhi 
with galaxy properties such as stellar mass, average stellar 
surface density and average gas surface density. We then elu- 
cidate the physical processes that determine the slope and 
scatter of these relations. Finally, we compare the predic- 
tions of the models with results derived from recent data 
sets. 

5.1 Model results 

In Fig. [7] we plot Mh 2 /M H i, Mh 2 /M*, Mhi/M* as a func- 
tion of stellar mass M, , mean stellar surface density /u* and 
mean gas surface density ^i gas for our model galaxies. In this 
analysis, we select the model galaxies with bulge fraction 
B/T < 15% in haloes with « vir >120km/s at z=0. These 
scaling relations should thus be regarded as appropriate for 
massive late-type galaxies. We do not consider early-type 
galaxies or the dwarf galaxy population in this paper, be- 
cause we do not think our disk formation model is likely 
to provide an accurate representation of how such galaxies 
assemble. 



We define fi, as 



where rso is the radius enclosing half the stellar mass of 
the galaxy. Note that we have not implemented a detailed 
model for how the stars in the bulge are distributed as a 
function of radius. In order to calculate rso in the models, we 
simply assume that the stellar mass in the bulge is included 
within r 5 o . ^ gas is denned as the mass- weighted mean surface 
density of the gas in the disk of the galaxy: 

Mgas = m gas /7Tfg as , (37) 

where the mass weighted mean radius of the gas disk f gaa is 
defined as 

/ rrfm gas 1 \ - , . , . 

rgas = Td^ — = ^—1*. r '( m 8»)i- ( 38 ) 

The sum in Equation (|38|) extends over the 30 rings used to 
represent the disk in our models. 

The left and right panels in Fig. [7] present results for 
the two H2 fraction prescriptions. The gray dots represent 
results for individual galaxies, while the blue curves with 
error bars indicate the mean value of the gas fraction (or 
molecular-to-atomic ratio) as a function of the scale pa- 
rameter on the x-axis. Error bars indicate the ±1<t scat- 
ter about the mean. The atomic gas fraction Mhi/M* is 
inversely correlated with stellar mass, stellar surface mass 
density and gas mass surface density. The molecular gas 
fraction Mh 2 /M 4 is correlated with the gas mass surface 
density, but shows little dependence on stellar mass or stel- 
lar mass surface density. The molecular-to-atomic gas ratio 
Mh 2 /Mhi is tightly correlated with the surface mass density 
of the gas, and more weakly correlated with stellar mass and 
stellar mass surface density. 

A comparison of the two panels in Fig. [7] shows that the 
two H2 fraction prescriptions produce only slightly different 
results. The most obvious difference is that the correlation 
between Mh 2 /Mhi and ^i gas is considerably tighter for the 
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Figure 7. Mhi/M*, Mh 2 /M, and Mg, /Mjji arc plotted as a function of mass weighted mean gas surface density fi gaa , mean stellar 
surface density (j,*, and stellar mass M* for model galaxies with B/T ^ 15% in dark matter haloes with i> v i r >120km/s at z=0. Results 
are shown for both H2 fraction prescriptions. In each plot, the gray dots are the model galaxies and the blue curves show the mean values 
of Mhi/M», Me 2 /M t , Mn 2 /Mm as a function of the scale parameters on the x-axis. The error bars indicate the ±lcr scatter about the 
mean and the bins have been chosen to contain an equal number of galaxies. We have also labelled each panel with the values of the 
mean scatter about the relation in logarithmic units (the scatter towards higher values is denoted "+", and the scatter towards lower 
values is denoted "-"). 



pressure-based prescription than for the prescription based 
on the Krumholz et al. (2009) models. The very tight re- 
lation between Mn 2 /Mni and ^i gas in the right hand panel 
arises because E gas dominates the contribution to the pres- 
sure for the vast majority of the late-type galaxies in our 
simulation. However, galaxies of given gas surface density do 
have a spread in metallicity, and this causes the increased 
scatter in the Mh 2 /Mhi versus fi sas relation in the left panel. 



5.2 Physical origin of the scaling relations 

We will now attempt to elucidate the physical origin of the 
scaling relations discussed in the previous section. Because 
the differences between the two H2 fraction prescriptions are 
small, we will only focus on the results from prescription 1 
in the following analysis. 

In our galaxy formation models, the mass and structure 
of the disk are set by parameters that cannot, in general, be 
directly observed. These are: 

(i) The virial mass of the dark matter halo M v i r . 

(ii) The halo spin parameter A (Equation 1 101) . 

(iii) The fraction of gas that has been accreted at recent 

epochs, /new 

In this paper, we define / ncw as the ratio of the mass of the 
gas that was accreted in the last Gigayear to the current 
total disk (stellar+gas) mass of the galaxy disk: 



./new — 



Wcooi(< lGyr) 



(39) 



m d i s k(z = 0) 

In Fig. [8] we plot the same scaling relations as in the 
left panel of Fig. [Jj but this time we divide each plane into a 



set of cells, and we bin galaxies according to which cell they 
occupy. We then plot the average value of M v i r , A and /new 
in each grid cell that is occupied by more than 100 galaxies. 
The results can be summarized as follows: 

(i) M vir : Panel (a) of Fig. |8] shows that more massive 
haloes host more massive galaxies with higher stellar surface 
mass densities and lower gas fractions. The reason is that 
more baryons are available to form stars in more massive 
haloes. In addition, feedback processes are less efficient at 
expelling baryons from the galaxy in massive haloes, so the 
fraction of the baryons that cool and turn into stars is also 
higher. 

(ii) A: Panel (b) of|8]shows that galaxies in haloes with 
larger spin parameters have higher values of Mhi/M* and 
Mh 2 /M«, but lower values of Mh 2 /Mhi- The spin parameter 
sets the contraction factor of the infalling gas. If it is large, 
a large fraction of the recently accreted gas will be located 
in the extended, low density regions of the outer disk, where 
a smaller fraction of the cold gas will be in molecular form. 
The resulting star formation rate surface densities will thus 
be low. This is why galaxies in haloes with large spin pa- 
rameter have low stellar surface densities, high gas fractions, 
but lower-than-average ratios of molecular-to-atomic gas. 

(iii) /new: Panel (c) of Fig. [8] shows that galaxies that 
have recently accreted a significant amount of new gas have 
higher values of Mhi/M* and Mn 2 /M t , but lower values of 
Mh 2 /A/hi- Because galaxies form from the inside-out in our 
model, gas that has been recently accreted will be located in 
the outer regions of the disks and have lower surface densi- 
ties and lower molecular fractions. To first order, therefore, 
the gas properties of galaxy that has experienced recent ac- 
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Figure 8. The relations between Mjji/M*, Mu 2 /M„, Mh 2 /Mhi vs. mass weighted mean gas surface density /i ga s, mean stellar surface 
density fit , and stellar mass M* from model results for disk galaxies at redshift=0, which reside in haloes with virial velocities greater 
than 120 km/s. In each grid, the colours represent the mean values of halo virial mass M vlr (subfigure a), spin parameter A (subfigurc 
b) and new gas fraction / ne w (subfigure c). Each grid cell contains at least 100 galaxies. 



cretion will be similar to a galaxy located in a halo with 
somewhat higher spin parameter. 

The degeneracy between A and /new can be analyzed 
in more detail by comparing panels (b) and (c). The basic 
trends along the y-axis are quite similar for both panels: 
galaxies with higher-than-average gas fractions and lower- 
than average molecular-to-atomic gas ratios could either 
have experienced a recent accretion event, or could simply 
have higher spin parameters. Some notable differences do 
emerge, however, in the plots of gas fractions and molecular- 
to-atomic gas ratio as a function of /% aB . Galaxies with 
Mgas > lO 8 M0kpc -2 with higher-than-average gas fractions 



are very likely to have experienced a recent gas accretion 
event. If one wishes to maximize the likelihood of identi- 
fying a galaxy that has experienced recent gas accretion, 
Fig. E] suggests that one should select: a) galaxies with high 
gas-to-star ratios and high gas surface densities, b) galaxies 
with high gas-to-star ratios and low molecular-to-atomic gas 
ratios. For case (b), one would need to seek additional evi- 
dence that the gas in the outer disk was accreted recently. 
We will discuss some ideas about how this might be done in 
the final discussion. 

The results in Fig. [8] can also help us understand the 
origin of the relations in Fig. (7) As seen in Fig. EI a), more 
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massive galaxies form in more massive haloes. Because feed- 
back effects are less effective at preventing the baryons in 
high mass haloes from cooling, a larger fraction of the avail- 
able baryons is converted into stars, leading to galaxies with 
higher stellar surface densities and lower gas-to-star ratios. 
This explains why Mhi/M* and Mh 2 /M, are anti-correlated 
both with stellar mass and with stellar mass surface den- 
sity. In Fig. _2 C ), we see that recent gas accretion rates are 
higher for low mass galaxies. This is again because supernova 
feedback prevents baryons from cooling effectively in low 
mass haloes. Higher recent gas infall rates lead to disks with 
lower molecular gas fractions. This explains why Me 2 /Mhi 
is positively correlated with stellar mass and stellar mass 
surface density. Because the trends in total gas-to-star ratio 
and molecular-to-atomic ratio work in opposite directions, 
there are tighter correlations between the atomic gas fraction 
Mhi/M* and stellar mass and surface density, than between 
the molecular gas fraction Mh 2 /M 4 and stellar mass and 
surface density. 

Fig. [7| shows that our models predict a very tight corre- 
lation between the atomic-to-molecular ratio and the mean 
gas surface density of the disk. This is simply a reflection 
of our H2 fraction prescriptions - the local surface density 
of the gas is the primary parameter that determines the 
fraction of gas that is converted to molecular form. Because 
galaxies with higher values of fj, sas have higher total gas frac- 
tions and higher molecular-to-atomic gas ratios, there is a 
clear positive correlation between Mh 2 /M, and Atgas- How- 
ever, because more of the gas is transformed into molecular 
form at high total gas surface densities, the relation between 
Mhi/M* and ^i gas is rather flat. 

5.3 Comparison with observations 

In this section, we will compare the disk galaxy scaling re- 
lations predicted by our models with observations. One of 
the major problems that we face is the lack of suitable data 
sets. Older data, for example as presented in the recent com- 
pilation by Obreschkow & Rawlings (2009), are beset by a 
variety of uncertainties and biases. The observed galaxies 
were selected using a variety of different criteria, and single 
dish observations did not often cover the entire disk of the 
galaxy. These data sets are therefore heterogeneous and the 
H2 masses are likely to be inaccurate; it is therefore unclear 
whether the scaling relations derived from them are useful 
for our purposes. For example, Figure 5 of Obreschkow & 
Rawlings (2009) appears to indicate that local galaxies have 
values of log(MH 2 /MHi) that are distributed relatively uni- 
formly from -f.5 to 1, i.e. if one believes these numbers, one 
would need to understand why galaxies exhibit a range in 
molecular-to-atomic gas ratio of more than a factor of 100! 
In contrast, our models predict a range in Mh 2 /Mhi closer 
to a factor of 10 (Fig. [7]). 

We turn once again to the small, but internally con- 
sistent THINGS/HERACLES data set. From the tables in 
Leroy et al. (2008), we derive and //gas from the stel- 
lar and gas surface density profiles according to the defi- 
nitions in Equation (|36|) & (|37[) . The comparison between 
data and models presented in Fig. [9] is very encouraging. 
The gas fractions and atomic-to-molecular gas ratios of the 
THINGS /HERACLES galaxies span the same range of val- 
ues as our model galaxies. In addition, one sees that: (1) The 
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relation between Mhi/M* and /x* and m* is steeper than 
that between Mh 2 /M* and /i* and m*, (2) The relation be- 
tween Mh 2 /M* vs /%as is tighter than that between Mhi/M, 
and ^igas, (3) Mh 2 /Mhi and (j, sas exhibit the strongest and 
tightest correlations. These results are all consistent with 
our model predictions. 

The current observations are not able to distinguish be- 
tween the two H2 fraction prescriptions. Adding more data 
points to Fig. [9] may not help very much, because the main 
effect on the scaling relations is a change in the scatter in 
some of the plots. One would have to understand the obser- 
vational errors in considerable detail to know which part of 
the scatter was real. Instead, we propose to find systematic 
effects in the scaling relations that arise as a consequence of 
the H2 prescription. The Krumholz et al. model predicts that 
much of the scatter arises as a result of metallicity differ- 
ences between different galaxies. This is illustrated in detail 
in Fig. 1101 where we plot the molecular-to-atomic fraction 
Mh 2 /Mhi as a function of /i gas and fj,*, and we colour-code 
each galaxy according to its gas-phase metallicity. As can be 
seen, H2 prescription 1 (the Krumholz et al. model) predicts 
that there should be a clear stratification in metallicity at 
fixed /igas (and to a latter extent at fixed p,), with the most 
metal rich galaxies having higher values of A/h 2 /Mhi. This 
is not seen for the pressure-based H2 prescription. Since the 
gas-phase metallicity can be estimated using strong emis- 
sion lines in optical spectra (e.g. Tremonti et al. 2004), this 
should be easily testable with future data sets, such as that 
provided by the CO Legacy database for GASS survey (see 
http://www.mpa-garching.mpg.de/COLD_GASS7l. 



6 SUMMARY AND DISCUSSION 

In this paper, we extend existing semi-analytic models to 
follow atomic and molecular gas in galaxies. We study how 
the condensed baryons in present-day disk galaxies are par- 
titioned between stars, HI and H2 as a function of radius 
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Figure 9. Model relations between Mh,/Mhi, Mh 2 /M„, Mjji/M* and stellar mass M*, mean stellar surface density /i* and mass 
weighted mean gas surface density /igas are compared with available observational data. The green dots are model results for disk 
galaxies at z=0 in dark matter haloes with v v i v > 120 km/s. The red circles have been derived from the observational results tabulated 
in Leroy ct al. (2008). 



within the disk. Our new model is implemented in the L- 
Galaxies semi-analytic code and is a modification of the 
models of Croton et al. (2006) and De Lucia & Blaizot 
(2007) , in which dark matter halo merger trees derived from 
the Millennium Simulation form the "skeleton", on which 
we graft simplified, but physically motivated, treatment of 
baryonic processes such as cooling, star formation, super- 
nova feedback, and chemical enrichment of the stars and 
gas. We fit these models to the radial surface density pro- 
files of stars, HI and H2 derived from recent surveys of gas in 
nearby galaxies, making use of data from SINGS, THINGS, 
HERACLES and the BIMA SONG surveys. 

We have used our models to explore how the relative 
mass fractions of atomic gas, molecular gas and stars are 
expected to vary as a function of global galaxy scale param- 
eters, including stellar mass, mean stellar surface density, 
and mean gas surface density. We have attempted to eluci- 
date how the trends can be understood in terms of the three 
variables that determine the partition of baryons in disks: 
a) the mass of the dark matter halo, which determines the 
total mass of baryons that is able to cool and assemble in 
the disk, b) the spin parameter of the halo, which sets the 
contraction factor of the gas, and thereby its surface density 
and molecular fraction, c) the amount of gas that has been 
recently accreted from the external environment. 

The main changes we have made to earlier models are 
the following: 

(i) Each galactic disk is represented by a series of con- 
centric rings. We assume that surface density profile of in- 
falling gas in a dark matter halo is exponential, with scale 
radius that is proportional to the virial radius of the 
halo times the spin parameter of the halo. As the Universe 
evolves, the dark matter halo grows in mass through merg- 
ers and accretion and the scale radius of the infalling gas 
increases. Disk galaxies thus form from the inside out in our 



models. The ring representation allows us to track the sur- 
face density profiles of the stars and g function of 
cosmic time. 

(ii) We include simple prescriptions for molecular gas 
formation processes in our models. We adopt two differ- 
ent "recipes": one based on the analytic calculations by 
Krumholz et al. (2008), in which /h 2 is a function of the 
local surface density and metallicity of the cold gas, and the 
other motivated by the work of Elmegreen (1989 & 1993), 
Blitz & Rosolowsky (2006), and Obreschkow et al. (2009), 
in which the H2 fraction is determined by the pressure of 
the ISM. 

(iii) Motivated by the observational results of Leroy et 
al. (2008), we adopt a star formation law in which Esfr oc 
Eh 2 in the regime where the molecular gas dominates the 
total gas surface density, and Esfr <x Eg as where atomic 
hydrogen dominates. 

Our work leads to the following conclusions: 

(i) A simple star formation law in which Esfr oc Eh 2 
leads to gas consumption time-scales in the inner disk that 
are too short. In this paper, we simply patch over this prob- 
lem by decreasing the efficiency of supernova feedback in the 
inner disk. 

(ii) The mean stellar, HI and H2 surface density profiles 
of the disk galaxies in our model are only weakly sensitive 
to the adopted H2 fraction prescription. The reason for this 
is that for typical L» disk galaxies, the local gas surface den- 
sity is the main controlling parameter for both recipes. At 
low gas surface densities, the H2 fraction depends sensitively 
on metallicity for the Krumholz et al. prescription, but con- 
siderably less sensitively on stellar surface density /i» for the 
pressure-based prescription. As a result, the correlation be- 
tween molecular-to-atomic fraction and fj, sas for local disk 
galaxies exhibits more scatter if the Krumholz et al. model 
is correct. 
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(iii) Our results indicate that galaxies that have re- 
cently accreted a significant amount of gas from the external 
environment are characterized by higher-than-average total 
cold gas content. If the galaxy has high gas surface den- 
sity, then this excess gas is an unambiguous signature of 
a recent accretion event, because the time-scale over which 
gas is consumed into stars is short in such systems. On the 
other hand, if the galaxy has low surface density, a higher- 
than-average total cold gas content could indicate a recent 
accretion event, but it may also mean that the galaxy has 
a higher-than-average spin parameter. Higher spin parame- 
ters result in disk galaxies with more extended distributions 
of cold gas, lower-than-average molecular-to-atomic ratios, 
and low star formation efficiencies. For these ambiguous sys- 
tems, one must seek additional evidence that the outer disks 
were assembled recently. 

Although these conclusions are somewhat open-ended, 
they do suggest avenues for further research. We believe that 
a more realistic way forward to solving the gas consumption 
timescale problem would be to model radial inflow of the gas. 
Attempts have been made to construct phenomenological 
models that do include radial mixing of the stars and gas 
in disks as well as the effect of this mixing on the chemical 
evolution of the stars formed in the solar neighbourhood (e.g 
Schonrich & Binney 2009). Results from hydrodynamical 
simulations also indicate that the gas tends to flow inwards, 
while the stars migrate outwards (e.g Roskar et al. 2008). 
The main way to distinguish between different scenarios may 
be the predicted metallicity gradients. We intend to explore 
these issues in more detail in future work. 

In our model results, although the surface density pro- 
files from the two H2 fraction prescriptions are very sim- 
ilar, the models indicate that one should, in principle, be 
able to confirm the metallicity-dependence of the molecular 
gas fraction predicted by the Krumholz prescription, if one 
measures the average gas-phase metallicities of nearby disk 
galaxies using emission lines. Alternately, one can break the 
degeneracy by observing systems where the metallicity is 
low but the pressure is high (Fumagalli, Krumholz & Hunt 
2010). 

Another interesting issue is whether a galaxy's location 
in the gas scaling relation diagrams can serve as a diagnostic 
as to whether it has accreted gas from the external environ- 
ment. Although the theory of gas accretion in galaxies has 
received considerable attention of late (e.g. Keres et al. 2005; 
Dekel & Birnboim 2006; Dekel et al. 2009), there is little di- 
rect observational evidence that this occurs in practice. This 
is true both for galaxies in the local Universe and at high 
redshifts, where gas accretion rates are expected to be much 
higher. Although average gas accretion rates are expected 
to be low at the present day, precise quantification of the 
expected scaling relations for equilibrium disk galaxies may 
allow us to identify a subset of systems which deviate sig- 
nificantly from the mean in terms of their gas content. Fol- 
lowing the conclusion (iii) , one may try to gain a better un- 
derstanding of the observationally detectable signatures of 
a recent gas accretion episode. Possible ways forward would 
be to look for signatures of recent accretion in the observed 
age gradients of the stars or in the metallicity gradients of 
the gas in the disk. One could also look for accretion sig- 
natures in the kinematics of the stars and the gas in the 
outer disks. Alternatively, one could search for evidence of 



complex structure (e.g. tidal streams or shells) in the stellar 
haloes of gas- rich galaxies (Cooper at al 2010). We intend 
to explore these possibilities in more detail in future work. 

Ongoing and future surveys, such as the Galex Arecibo 
Sloan Survey (GASS) (Catinella et al. 2010) and the COLD 
GASS survey carried out at the IRAM 30m telescope (Sain- 
tonge et al. in preparation) will enable us to quantify the 
scaling relations discussed in this paper in considerable de- 
tail. These surveys will provide interesting targets for follow- 
up programs, which may help us understand that extent to 
which galaxies still accrete gas at the present day. In the next 
few years, it will become possible to observe gas in galaxies 
at higher redshifts using facilities such as ALMA and Square 
Kilometer Array pathfinder experiments such as ASKAP or 
MEERKAT. We are certain that our simplified treatment 
of disk formation in concentric rings that undergo no radial 
mixing will not be a good way to describe the assembly of 
the clumpy, highly turbulent disks that are now known to 
exist at z ~ 2 (e.g. Genzel et al. 2008). Nevertheless, we be- 
lieve that our models may still be useful in elucidating the 
gaseous and chemical evolution of disks over a somewhat 
smaller range in lookback time. 
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